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TRANSFER- FUNCTION- PARAMETER ESTIMATION FROM FREQUENCY 
RESPONSE DATA - A FORTRAN PROGRAM 
by Robert C. Seidel 
Lewis Research Center 

SUMMARY 

A FORTRAN computer program designed to fit a linear transfer function model to 
given frequency response magnitude and phase data is presented. A conjugate gradient 
search is used that minimizes the integral of the absolute value of the error squared 
between the model and the data. The search is constrained to insure model stability. 
Scaling of the model parameters by their own magnitude aids search convergence. Ef- 
ficient computer algorithms are used for the calculation of the cost function and gradi- 
ent and execution of the search. This results in rapid calculations and a small pro- 
gram size suitable for running on a minicomputer. Search convergence times for 
different model structures and parameters estimates for a sample problem derived 
from experimental data are presented. 


INTRODUCTION 

Frequency response techniques have long been a fundamental tool of system design 
and testing. One problem often faced is that of fitting a transfer function to given am- 
plitude and phase data. One of the most commonly used methods for doing this is the 
Bode plot graphical technique (ref. 1). Analytical techniques (refs. 2 to 4) that fit 
rational polynomial transfer functions using the digital computer have also been used. 

In reference 5 the mathematical framework is presented for a gradient search solution 
which, not being limited to rational polynomials, can handle dead time parameters, but 
there are no computer programs or practical results presented. 

Reported herein is a FORTRAN computer program designed to fit a transfer func- 
tion to frequency response data. The program has been used to obtain the transfer 
function models reported in references 6 and 7 for supersonic mixed- compression inlets. 
Search convergence times for different model structures and parameter estimates for a 



sample problem derived from reference 6 data are presented. A scaling rule is incor- 
porated to improve search convergence, and a constraint is placed in the search to in- 
sure model stability during the search. Efficient algorithms used for several different 
calculations result in a faster program with reduced storage requirements. 

The program formulation is first presented, including formulation of the cost func- 
tion, the model transfer function, the gradient calculation, and the techniques for scal- 
ing and constraining the model parameters. The computer program is discussed, and 
a sample problem illustration is given. In appendices A to D the report symbols, the 
computer program, the sample problem computer listing, and the proofs for various 
algorithms are presented, respectively. 


GENERAL PROGRAM FORMULATION 
Cost Function 

The desired result is to match a model transfer function to given plant magnitude 
A(w) and phase 9 ( w) frequency response data over a range of frequencies w. Ex- 
pressed in the form of a complex number A(w)e^ w \ the plant frequency response 
should be equivalent to the model transfer function G(jw,_b) where jw is the Fourier 
variable and b is a vector of model parameters. A cost function J(b) representing the 
fit error between the model and the plant data is formulated to be the integral of the 
absolute value of the error squared between the model and the data. 

\ 

J(b) = J | G(jw, b) - A(w)e^ w ^| dw (1) 

In the computer pregram all calculations are in the frequency domain. However, 
as an interpretation of equation (1), Parseval's theorem can be used to transform it 
into the time domain. As shown in reference 5, the result is that the cost function is 
the integral of the error squared between the plant and model impulse responses. 

The equation (1) is a continuous integral and is not suitable for evaluation with a 
digital computer; it is also not suitable for tabulated frequency response data. If the 
tabulated data are a sufficiently representative sample of the plant, the integral can be 
approximated by replacing the integration with a piecewise linear summation (trapezoi- 
dal rule): 
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1 
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je(w.) 

G(jw i ,b) - A(w i )e 


(w i + l - w i-l> 


( 2 ) 


where Wq - w^, + 1 = Wjj , and > 2 is the number of data points in the fre- 

quency response. Equation (2) is not the usual form of the trapezoidal rule. It has been 
modified to increase the program speed. The steps in obtaining the equation (2) form 
from the trapezoidal rule are presented in appendix D. 

The choice of the frequencies for the tabulated data can be important. The fre- 
quency interval (w i+1 - w i _ 1 ) weights the fit error at each frequency of equation (2). If 
the frequencies are widely spaced, the value of the error at a particular frequency could 
have a disproportionate effect. In certain cases it might be desired to favor a close fit 
between the model and the plant over a portion of the frequency range. The error could 
be purposely weighted to assign a relative importance to the fit as a function of fre- 
quency. This could be done to reflect uncertainties in the data, for example. The com- 
putation time increases proportionately with N d< Therefore, the user may wish to 
select a limited number of data points to represent the plant. This is the likely case 
for pulse testing followed by a fast Fourier transform of the data, which produces a 
large amount of closely spaced frequency response data. The measurement noise on 
actual systems has not been too favorable for pulse testing to date, however. The single 
frequency sinusoidal testing, by comparison, concentrates power at each frequency, 
virtually eliminating serious measurement noise problems. 


Model Transfer Function 

The model transfer function is assumed to have a linear structure expressed in 
factored form. This is the general form familiar to Bode plot users. In this form the 
transfer function G(s,b), where s is the Laplace variable, consists of the product of 
a number of possibly different parameter factors. 


G(s, b) = s 


m -n 
k P <1 
K 1 


1=1 


g m (s >b) £ m = 1,2, 3, 4, 6, 8 


(3a) 


The g m (s,b) of equation (3a) are the following parameter factors: 
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(3b) 


g (s,b) = — + 1 
‘ bl k 


g 2 (s,b) = b2 


(3c) 


g 3 ( s ,b) =( + l 1 

\ b3 k 


v-1 


s 2 2sb4 k 

g 4 (s,b) = — + 

b5 k b5 k 


+ 1 


-1 


g 6 (s,b) = 




.b7 2 b7 k 


(3d) 


(3e) 


(3f) 


g 8 (s,b) = e -sb8 (3g) 

The factors of equations (3b), (3d), (3e), and (3f) are assumed to be repeatable in 
equation (3a). The elements of the b vector are the zeros bl k , the gain b2, the poles 
b3 k , the quadratic zero damping ratios b4 k , the quadratic zero natural frequencies 
b5 k , the quadratic pole damping ratios bf^, the quadratic pole natural frequencies b7 k , 
and the exponential (dead time) b8. The number of g m (s,b) factors in the transfer 
function is the number of parameters in the b vector m p minus the number of quad- 
ratic factors n^. The exponent of the free s is kj. For example, a particular 
G(s,b), where 


b T = (b2,b8,bl 1 ,bl 2 ,b6 1 ,b7 1 ) 


is 


G(s,b) 



(4) 
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Gradient Calculation 


The gradient of J(b) with respect to b is information required by a gradient search 
to find parameters b that minimize J(b). The gradient can be expressed analytically 
by differentiating the equation (2) cost function. The result can be expressed in a form 
that facilitates rapid computation by invoking the following identiy for a complex func- 
tion E: 


(|e| 2 )' = (EE*) * = E(E ')* + E 'E* =2/?e (E(E')*) (5) 

Here, the prime indicates a derivative and the asterisk indicates a complex conjugate. 

A proof of this identity is presented in appendix D. The gradient of J(b) is 


Nj 


VJ(b) =/?e |G(jw.,b) - A(w i )e 

U1 


je(w.)' 


VG(jw.,b)*(w. +1 - w^j) 


( 6 ) 


Furthermore, the VG(jw., b) factor can be written to include all the possible combina- 
tions of factors for G(jWj, b) by defining a vector Z(s, b) such that 

VG(s,b) = G(s,b)Z(s,b) (7) 

By using G(s,b), which is already computed, the expression is computationally effi- 
cient. The Z(s,b) elements are internal to the computer program which assigns to 
each parameter the proper form. The details of the Z(s,b) forms are presented in 
appendix D. 


Parameter Scaling 

It was soon apparent from initial trials without parameter scaling that many 
searches would not converge. Gradient search techniques generally require parameter 
scaling to obtain efficient search convergence. An ideal scaling rule would obtain nearly 
spherical cost functions (where the cost is more or less equally sensitive to each 
parameter of b) and would tend to result in rapidly converging searches. The scaling 
law adopted applies to the specific form of the factored transfer function parameters. A 
rule of thumb scaling law was incorporated such that parameters tend to have equal in- 
fluence. This law scales each parameter by its own magnitude. To show this, a new 
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scaled parameter vector p with elements p^ is defined 


p. = — - — b. 

' H ‘ 


( 8 ) 


where the (b^ is the magnitude of the b vector elements b i at the start of each line 
search. The elements of the scaled gradient J(p) are formed using the chain rule and 
equation (8): 


_§J = ^i _3J = / b \ _3J 
3Pj 8p i 3b^ i h 3^ 

The gradient is computed only at the start of each line search. At this time b. = (b.) . 
Thus, the scaled gradient terms are merely the unsealed gradient terms multiplied by 
their corresponding parameter values. 

The p vector does not actually appear in the computer program. To avoid a con- 
tinual multiplication of the scaling factor times its scaled parameter and gradient, the 
search subroutine was specialized to handle the actual parameters and the scaled gra- 
dient. 

Normally, a conjugate gradient search is rescaled only after completing a number 
of line search iterations equal to at least the number of parameters plus one. But re- 
scaling here is done after each iteration and is a nonstandard modification to the con- 
jugate algorithm of reference 8. The advantages for the search program are less logic 
and storage requirements for scaling parameters to their current value. In tests there 
were no significant differences found in either the number of iterations or convergence 
accuracy between rescaling every iteration and after the usual number. This is be- 
lieved to result from the fact that away from the minimum J(b) is likely to be a poor 
approximation to the theoretically required quadratic b costs and near the minimum 
scaling changes are small. 


Parameter Constraints 

Ideally, it would not be necessary to constrain the model parameters for stability 
as the infinite error cost function boundary separating stable and unstable models would 
be reflected in the cost function J(b). However, large changes in b during the search 
may enable a jump across a relatively narrow error boundary. The solution to this 
problem was to force a sharp step size reduction on parameters trying to cross zero. 
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Zero is critical because it is the point of instability for factored parameters. A modifi- 
cation to the gradient search algorithm constrains the parameters to remain the same 
sign during the search. This is done by assigning to a parameter about to cross zero a 
small value on the same side of zero. Thus, an initially stable model must remain 
stable. This action also prevents changes between right and left hand plane zeroes, and 
it also prevents changes between dead and lead time exponentials. 


COMPUTER PROGRAM AND SAMPLE PROBLEM 

The FORTRAN IV computer program reported here mechanizes the transfer func- 
tion parameter identification. The usual procedure to identify experimental plant data 
consists of entering tabulated plant amplitude A(wp and phase shift 0(w.) frequency 
response data and a trial factored form transfer function. The user also supplies 
roughly estimated parameter values and a code number to identify them as time con- 
stants, gains, etc. The program can then be directed to iterate on the supplied param- 
eter values until it has optimized the fit to the data. The program output contains the 
optimized parameter values and a number indicative of the fit error J. Additional out- 
put includes the model frequency response. The computer program will not change the 
signs of the user supplied parameter values nor will it change the given transfer function 
structure. 

The computer program was written primarily for a conversational (time -sharing) 
computer, but it could also be run on a batch process computer. The effort to keep 
down program size and increase program speed was to enable its use on a minicomputer. 
The basic program was run on an SEL 810-B minicomputer in a nearly on-line identifi- 
cation implementation. The intent to obtain a small program package has ruled against 
writing a program with such conveniences as a plotting package, error detection rou- 
tines, and alternative search routines. The conjugate gradient search was chosen here 
because of its simplicity and ease in rescaling and constraining the parameters. Other 
search routines may converge in fewer iterations. But this could be overbalanced by 
the resulting increase in computing complexity and scaling and constraining problems. 

It has not been required to test the program on higher order models (greater than 
15 parameters), but generally long convergence times are indicated. 

There are provisions in the computer program to make it more general. On occa- 
sions it may be desired to not iterate on known parameters in the model, but to leave 
them fixed. This can be accomplished in the program through the computer variable 
JD, which is the number of fixed variables. Another option is helpful for reducing 
higher order analytical models to lower order models. Known plant transfer functions 
can be entered (rather than tabulated plant frequency response data) by setting the com- 
puter variable INST to 5. 
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The computer program may be used to design a series compensator that will make 
a plant look like any desired transfer function. For example, the desired transfer func- 
tion for the plant plus compensator divided by the plant may be entered or built. Then 
the compensator assumes the role of the model transfer function. 

The program package consists of a main program called MODEL and two subrou- 
tines called CGFM and CFG. The MODEL program organizes the data and is written 
primarily for a time- sharing type of computer. The subroutine CGFM conducts the 
conjugate gradient search. The subroutine CFG does the actual work of computing the 
cost function and gradient. These programs are described in detail in appendix B. 

The problem is to fit transfer function models to given experimental frequency re- 
sponse data using the computer program. The data are the response of an inlet- 
pressure signal to a bypass-door control area and is taken from reference 6. The data 
are the points plotted in figure 1. In handling experimental data the user usually does 
not know beforehand the model structure. The most useful analytical model may be a 
compromise between lower order and accuracy of fit to the data. In trying to fit the 
data a user may iteratively work up to more complex models. The dashed curve of fig- 
ure 1 is the fit of the following three- parameter model to the data: 


G(s, b) = 


|_17. 5(2)7) 


■+ 1 


t 2s(0. 262) ] j 
L130(2 tt)J 130(2tt) 


(10) 


The solid curve in figure 1 is the fit of a six-parameter model to the data. A computer 
listing of this solution is presented in appendix C. There it is shown that the initial 
estimate of 


G(s,b) 


— + 1 

.1(277) , 


s + 1 

S +1 

s +1 1J 

IT s 

2 . 2s w , A 

.1(277) . 

.10(2)7) 

1000(2)7) J] 

[[.1000(277 )_ 

1000(277) J 


(ID 


converged to the optimized model 


G(s,.b) 


2 + 1 

(4. 91)(277 ) J 


S + ll 

s +1 

S i 1 J 

IT s 

2 , 2s(0. 265) , j] 

_(3.70)(2t7) 

(27. 0)(2tt) 

(309)(2)7) J] 

[L(145)(2jt)_ 

(145)(2t7) J 


( 12 ) 
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after 93 iterations. The six-parameter fit is visibly better over the frequency range of 
interest, but both are optimum for their given structure. 

The initial parameter estimates can affect the search convergence time and some- 
times result in detecting a local rather than a global minimum. The convergence times 
and number of iterations for initial parameter estimates that were off by about a factor 
of 10 are listed in table I for four different model structures. Convergence times 
ranged from 24 seconds for the five -parameter model to 1.8 seconds for the two- 
parameter model (CPU time IBM 360-67 TSS). The fit error J(b) as calculated by the 
program for various order models are tabulated in table I. The error ranges from 5. 13 
for the two-parameter model to 0. 59 for the six-parameter model. The program dis- 
carded zeroes for models with less than four poles by making them very large. Like- 
wise, the program discarded a four-pole model by making the fourth pole very large. 
Two simple poles converging to the same value indicated that a quadratic pole would 
probably do better. 


CONCLUSIONS 

A FORTRAN computer program is presented which fits a linear factored form 
transfer function to given frequency response data. The identification process is often 
an important first step in system control design. The computer program is based on a 
conjugate gradient search procedure that minimizes the error between the given 
frequency-response data and the frequency response of a transfer function that is sup- 
plied by the user. The user- supplied transfer function consists of a product of terms 
that can include a gain, an integrator (or differentiator), a dead time, multiple first- 
order leads and lags, and multiple second- order leads and lags. 

The usual procedure for use of the program consists of entering a table of amplitude 
and phase-shift frequency-response data and a trial, factored form transfer function; 
the user also supplies roughly estimated parameter values such as time constants, 
gains, etc. The computer program then iterates on the supplied parameter values until 
it has optimized the fit to the experimental data; the excellence of this fit is constrained 
by the form of the user- supplied transfer function. The computer output consists of the 
optimized parameter values and a number indicative of the fit error. The computer 
program does not change the signs of the user-supplied parameter values, nor does it 
change the given transfer- function structure. The user can "build" a transfer function 
by trying a series of increasingly complex transfer -function structures being guided by 
the fit errors in a systematic manner as illustrated by the sample problem. 

An alternative use is to reduce the complexity of a high-order transfer function that 
might result from a dynamic analysis. An option is provided in the program for enter- 
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ing the complex transfer function in factored form to take the place of the normal table 
of frequency response data. The simpler model can then be "built” by the usual pro- 
cedure until it is a reasonable fit to the more complex model. 

The conjugate gradient search routine was chosen for the error function minimiza- 
tion because of its simplicity and ease in constraining parameters. In the program, 
scaling the parameters by their own magnitude was demonstrated to result in good 
search convergence. Constraining the parameters to remain the same sign during the 
search insures model stability. The use of time and space saving algorithms resulted 
in a small and fast program. The program as presented is written primarily for a con- 
versational (time sharing) computer, but it could be run on a batch process computer. 
The basic program has been run on an SEL 810-B minicomputer as part of a nearly on- 
line transfer function identification implementation. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, June 23, 1975, 

505-05. 
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APPENDIX A 


SYMBOLS 

A system frequency response magnitude 

b parameter vector, nip x 1 

it. 

bl, b vector k Ln zero, rad/sec 

" k l 

b2 b vector gain parameter, (rad/sec) 

it. 

b3 k b vector k in pole, rad/sec 

IL 

b4^ b vector k n quadratic zero damping 
th 

b5 k b vector k in quadratic zero natural frequency, rad/sec 

1L 

— vector k quadratic pole damping 
b7k b vector k* k quadratic pole natural frequency, rad/sec 
b8 b vector exponential (dead time), sec 
c real variable 

d real variable 

E complex function, E = c + jd 
f real function 

f, f(Xj) 

G model transfer function 

g transfer function factored parameter factor 

J cost function 

j imaginary number, yTl 

jw Fourier transform variable, rad/sec 

kj integer, exponent of free s's in model 

mp integer, number of elements in b vector 

N d integer, number of frequency points in numerical integration 

n integer, number of quadratic factors 

VI 

p transformed parameter vector, m p x 1 

s Laplace variable, sec"* 
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w 

frequency, rad/sec 


w i 

frequency at i^ data point, rad/sec 

x i 

1L 

variable at i data point 


z 

partial product in gradient vector, 

m p xl 

9 

system frequency response phase, 

deg 

Superscripts: 


f 

derivative 


* 

complex conjugate 


T 

transpose 


Subscripts: 



value at start of iteration 


i 

i** 1 element 


£ 

jT n factored parameter factor 


m 

type of factored parameter factor 
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APPENDIX B 


COMPUTER PROGRAM 

The FORTRAN computer program package for transfer function parameter identifi- 
cation consists of the main program MODEL, subroutine CFG to calculate the cost func- 
tion and gradient, and subroutine CGFM to conduct the conjugate gradient search. Sym- 
bols used in the programs are defined separately for each program. The programs are 
dimensioned for a frequency point maximum of 50 and mp model parameter maxi- 
mum of 15. The vectors dimensioned for N^ are AMP, PHA, and PLANT. The 
vectors dimensioned for m are B, G, ID, and Z. The vector H is dimensioned 2m 

r r 

and the vector W is dimensioned N^ + 1. 


Program MODEL 

The program MODEL handles the data and converts system frequency response data 
to complex numbers. The computer variables in the program are defined in the section 
Program MODEL Variable List. A flow chart for the program MODEL is presented in 
figure 2. 

The program starts by printing a heading referencing the namelist variables and 
variable codes. Then the program prompts for NAM3 namelist data. The namelist 
variables are entered according to the rules for entering FORTRAN namelist data. The 
NAM3 variables are AMP, ND, PHA, W, and KPR. The program prints the NAM3 
variables if KPR = 1 and then prompts for NAM2 namelist data. The NAM2 variables 
are B, Kl, ID, LIMIT, N, and JD. Then the NAM2 variables are printed and a prompt 
for the INST variable is printed. 

The INST variable is entered in II format. The INST code values are 1, 2, 3, 4, 

5, and 6. Making INST = 1 causes a search for the optimum parameters. After the 
search results are printed, another INST prompt is issued. Making INST = 2 returns 
the program to request NAM2. Making INST = 3 returns the program to request NAM3. 
Making INST = 4 causes the model response to be printed. The frequency is printed 
under W, the model frequency response amplitude under AMP, and phase under PHA. 
Then the zero iteration search results are printed and another INST prompt is issued. 
Making INST = 5 is the same as making INST = 4 except that the PLANT complex num- 
bers are set equal to the model. This option is used to save entering AMP and PHA 
data and to build an analytical model for the system. Making INST = 6 stops the pro- 
gram. 

The frequencies W are input in radians, but they may also be input in hertz if the 
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B vector is appropriately scaled by 2n (see B vector in Program MODEL Variable 
list). 

The following is a FORTRAN listing for the program MODEL 1 : 


C PROGRAM MODEL FOR. TRANSFER FUNCTION IDENTIFICATION 

DIMENSION B( 15) , G( 15) , I D( 15 ) , AMP( 50 ) , PHA( 50 ) , W( 51),H(30) 
COMPLEX PLANT (50) 

EXTERNAL CFG 
LOGICAL KG 

COMMON/ I DTN/ PLANT, W, I D, Kl, I NST, ND, KNT, NP 
COMMON /FMC/KOUNT, KG 

NAMELIST /NAM3/W, AMP, PHA, ND, KPR /NAM2 /L I Ml T, N, Kl, B, I D, JD 
WR I TE( 6. 110) 

110 FORMAT ( ' NAM3=( W, AMP, PHA, ND, KPR) ' , / - 

1, ' NAM2 = ( B, I D, Kl, L I Ml T, N, JD) ' , / - 

2, ' ID=(1=Z;2=G;3®P;4,5=CZD,CZW;G,7=CPD,CPW;8=DT)',/- 

3, ' I NST®( 1=SEARCII, 2-HAM2, 3=NAM3, 4 s PR I NT, 5 = PLANT “MODEL, 6 = - 
4STOP)',/,' IER=(0=CONV, l=NOT CONV, 2-ERROR) ' ) 

10 WRITE (6,180) 

180 FORMAT (' NAM 3? ' ) 

READ (5, NAM3 ) 

IF(KPR.EQ.l) WR I TE(6, 310) 

310 . FORMAT ( 6H PLANT , /, 8X, 1UW, 14X, 3HAMP, 11X, 3HPHA) 

120 FORMAT (1P3E15.5) 

DO 20 J«1,ND 

IF(KPR.EQ.l) WR I TE( 6, 120 ) W( J ) , AMP( J ) , PHA( J ) 
TI1ETA=PHA(J)*. 0174 532 9 2 
20 PLANT ( J) =AMP( J )*CEXP(CMPLX(0., THETA) ) 

W( ND+ 1 ) »W( ND ) 

30 WRITE (6,181) 

181 FORMATC NAM2? 1 ) 

READ ( 5, NAM2 ) 

NP=N+JD 

WR I TE ( 6, 130 ) LIMIT, N,K1, JD, ( I D( J), J-l, NP) 

130 FORMATC L I Ml T, N, Kl, JD= ' , 4 I 3, ' I D= ' , 1 5 I 3 ) 

GO TO 87 

777 WRITE (6,182) 

182 FORMATC INST?') 

READ ( 5, 140) I NST 

140 FORMAT (II) 

GO TO (40, 30, 10, 50,45, 1000), INST 
45 WR I TE( 6, 145) 

145 FORMATC PLANT SET EQUAL TO MODEL') 

50 WR I TE(6, 160 ) 

160 FORMAT ( 6H MODEL, /, 8X, 1HW, 14X, 3HAMP, 11X, 3HPHA) 

40 SIZE®. 1 

EPS=1. E-5 

I TER=L I Ml T*( 1- I NST/4 ) 

KNT = 0 

CALL CGFM(CFG,N,B,F,G, SIZE, EPS, ITER, I ER, H ) 

WR I TE( 6, 150 ) F, I ER, KOIJNT, KNT, SIZE 

150 FORMAT ( ' F= ' , 1PE16 . 7, ' I ER= * , I 2 , ' KOUNT, KNT, S I ZE* ' , - 
1 13,14, 1PE10 . 3) 

87 WRITE (8,151) ( B( J) , J-l, NP) 

151 FORMATC B® ' , 1P8 Ell . 4) 


1 


A 


dash at end of a line 


signals continuation card follows. 
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1000 


GO TO 777 

STOP 

END 


AMP 

B 

CFG 

EPS 

F 

G 

H 

ID 


IER 

INST 


ITER 

J 

JD 

KG 

KNT 

KOUNT 


Program Model Variable List 

system frequency response magnitude A(w), vector (input variable) 

model parameter b, vector (input variable). Damping ratios must immedi- 
ately precede their natural frequency terms. Constant (JD) parameters 
must follow last. 

subroutine (declared external in model program) which computes the cost 
function and gradient 

parameter change defining search convergence, e.g. , 10“® 

F = J (b), cost function 

cost function (scaled) gradient, vector 

storage, vector 

integer vector that identifies corresponding parameter in B as to type (input 
variable): 1 = zero, 2 = gain, 3 = pole, 4 = complex zero damping, 5 = 
complex zero natural frequency, 6 = complex pole damping, 7 = complex 
pole natural frequency, 8 = dead time 

search convergence parameter: 0 = search convergence (SIZE reduced to 
EPS or less within LIMIT iterations), 1 = search not converged, 2 = prob- 
able error occurred 

branching instruction parameter (input variable): 1 = search for optimum, 

2 = return to NAM2 namelist, 3 = return to NAM3 namelist, 4 = print fre- 
quency and model frequency response amplitude and phase, 5 = set con- 
verted system PLANT frequency responses equal to model values and print, 
6 = stop 

set equal to LIMIT except for the INST = 4 or 5 in which case ITER = 0 
index of element in vector 

number of noniterative constant parameters in model (input variable) 
logical variable .TRUE, means compute gradient 
count of cost function evaluations 
count of line search iterations 
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KPR 

K1 

LIMIT 

N 

ND 

NP 

PHA 

PLANT 

SIZE 

THETA 

W 


if equal to 1 causes NAM3 variables to be printed (input variable) 

exponent of free s in model (input variable) 

maximum number of iterations (input variable) 

number of iterative parameters of B in model (input variable) 

number of frequency points over which integration is performed 

total number, iterative plus constant, of model parameters, NP = N + JD 

system frequency response phase 0(w), in deg vector (input variable) 

converted system response to equivalent complex numbers vector 

parameter step size; e. g. , set to 0. 1 at start of search 

scratch variable for radian phase 

radian frequency w, vector (input variable); can be input in hertz instead of 
radians but then the units of parameters in B may change. (The B 
parameters with ID numbers of 1, 3, 5, and 7 need to be divided by 2 n to 
obtain their corresponding value in hertz, a B with an ID = 8 needs to be 
multiplied by 2n, and B’s with ID = 2, 4, and 6 are not affected) 


Subroutine CFG (N, B, F, G) 

The purpose of this subroutine is to compute the cost function and gradient. The 
gradient G is the scaled gradient VJ(p). The unsealed gradient for other search pro- 
grams could be obtained by dividing each term in G by the corresponding parameter 
in B. 

The CFG variables are defined in the CFG Variable List. Those variables carried 
over in the common block and subroutine call are labeled the same as those described 
in the program MODEL and are not repeated again. Figure 3 is a flow chart of CFG. 


FORTRAN Listing for Subroutine CFG 


SUBROUTINE CFG(N,B,F,G) 

DIMENSION B(1),G(1), ID(15),W<51) 

COMPLEX PLANT(50),Z(15),S,ST,GM,E 
LOGICAL KG 

COMMON/ I DTN/PLANT,W, I D, Kl, I NST,ND, KNT,NP 
COMMON/ FMC/KOUNT, KG 
KNT-KNT+1 
DO 30 J=1,N 
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30 


G(,l)*0. 

WS-W(l) 

F =0 . 

DO 100 J= 1, ND 
S=CMPLX ( 0 . , W( J ) ) 

GM»S**K1 
DO 8 K=1,NP 
I DK= I D( K) 

GO TO (1,3,1, 2, 8, 2, 8,4), I DK 

1 GM=GM* (S/B(K) + 1.)**(2-IDK) 

I F ( KG ) Z(K)=S/(S+B(K))*( FLOAT ( I DK)-2. > 

GO TO 8 

2 ST=(S/B(K+1))**2+2.*S*B(K)/B(K+1)+1. 
GM=GM*ST**( 5-1 DK) 

I F( . NOT. KG) GO TO 8 

Z(K)=2.*S*B(K)/( B(K+1)*ST)*(5. -FLOAT ( I DK) ) 
Z(K+1)=-Z(K)*(1.+S/(B(K)*B(K+1))) 

GO TO 8 

3 GM=GM*B(K) 

Z(K)=CMPLX(1.,0.) 

GO TO 8 

4 GM=GM*CEXP(-S*B(K)) 

Z(K)=-S*B(K) 

8 CONTINUE 

E»GM-PLANT ( J ) 

F-F + (W( J+1)-WS)*REAL( E*CON.JG(E) ) 

I F( .NOT. KG) GO TO 100 
ST=(W( J+l)-WS)*E*CONJG(GM) 

DO 80 K=l, N 

80 G( K)=G( K) + REAL(ST*CONJG(Z( K) ) ) 

I F( INST.EQ. 1) GO TO 100 
AM=CABS( GM) 

DEG=ATAN2(AIMAG(GM) # REAL(GM))* 57. 295780 
WRITECB.IZO) W(J),AM,DEG 
120 FORMAT (1P3E15.5 ) 

I F( INST.EQ. 5) PLANT ( J ) =GM 
100 WS=W(J) 

RETURN 

END 


Subroutine CFG Variable List 

AM |G(jw.,b)| 

DEG ZG(jw.,b), deg 

E error between model and data 

GM partial product becomes G(s,b) 

IDK ID(K) value 

J frequency index 

K parameter index 



s 


s 


ST temporary value 

WS saved W value 

Z partial product in model gradient, vector 


Subroutine CGFM (CFG, N, B, F, G, SIZE, EPS, ITER, IER, H) 


The purpose of subroutine CGFM is to perform the conjugate gradient search func- 
tion minimization. Several nonstandard modifications relative to the conjugate gradient 
search described in reference 8 exist in the CGFM subroutine. The gradient is the 
scaled gradient VJ(p) and every iteration updated b parameters change the scaled co- 
ordinate system. Another nonstandard modification to the search is that the signs of 
the parameters b are not allowed to change during the search. 

The minimum for each line search is found by either doubling or taking one third of 
the parameter values. The process is continued until the minimum cost is known to be 
bracketed. Then the minimum point is obtained from the minimum of a quadratic curve 
fit to three points. The gradient is not used in this process. To increase the program 
speed, there is logic to signal whether the gradient is needed; if not, only the cost 
function value is required of the CFG subroutine. 

The CGFM program variables are defined in the CGFM program variable list. 
Variables carried over in the common block and subroutine call are labeled the same as 
those described in the main program and are not repeated again. Figure 4 is a flow 
chart of CGFM. 


FORTRAN Listing of Subroutine CGFM 


SUBROUTINE CGFM( CFG, N, B, F, G, S I Z E, EPS, ITER, I ER, H ) 
DIMENSION B(1),H(1),G(1) 

LOGICAL KG 
COMMON/ FMC/KOIJNT, KG 
C INITIALIZATIONS 
KOUNT=0 
STEP=2 . 

I ER=-5 
5 BETA=0 . 

NCYC-0 
15 K=0 

DO 20 J= 1, N 
20 H(J)=B(J) 

KG=.TRUE. 

CALL CFG( N, B, F, G) 
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C TEST FOR STOPPING SEARCH 

I F( KOUNT . GE. ITER) I ER=1 
I F(SIZE.LT.EPS) I ER=0 
40 I F( I ER. GT. -2 ) RETURN 
C COMPUTE PAST GRADIENT WEIGHTING 
TSQR=0 . 

DO 50 J = 1,N 
50 TSQR=TSQR+G( J)**2 

I F(NCYC.EQ.O) GO TO 60 
BETA=TSQR/TSAVE 
60 SCALE=0 . 

DO 70 J=»1,N 
JH= J + N 

H(JN)=-G(J)+BETA*H(JN) 

70 SCALE=SCALE+ABS(H( JN) ) 

I F( SCALE . GT . 0. ) GO TO 80 
I ER=0 
GO TO 40 

80 SCALE=SIZE/SCALE 
TSAVE=TSQR 
NCYC=NCYC+1 
FSS=F 
L=1 

C UPDATE B'S 

100 DO 110 J=1,N 
JN=J +N 

110 B( J)=H( J)*ABS( 1 . +SCALE*H( JN) ) 
FS=F 


KG=. FALSE. 


CALL CFG( N, B, F, G) 

C LOGIC^ CHANGES UNE^SEARCH STEP SIZE OR CONCLUDES SEARCH 

IF(F.LT.FS) GO TO 130 
I F(L.GT.l) GO TO 140 
GO TO 150 

120 IF(F.LT.FSS) GO TO 140 
GO TO 160 
130 SCASV = SCALE 
L= L+l 

SCALE= SCALE*STEP 
FSS=FS 

I F(L. LT.15) GO TO 100 
I ER=2 
GO TO 40 

i C ( n F ' T nH U f, D . Rl AT . 1 ? ( ( ' URVE T0 3 PTS * JACKETING LINE SEARCH MIN. 
140 DO 148 J»l / N 

JN=J +N 


147 

148 


R1=H( J) 

R2=H( J)*ABS( l.+SCASV*H( JN) ) 

R3=B( J) 

I F( L . GT. 3) R1=R1+(R2-R1) /STEP 
X1=(FSS-FS)*(R1-R3) 

X2=(FSS-F)*(R2-R1) 

I F(ABS(R2-R3).GT. EPS/4.) GO TO 147 
IF(L.GT.l) B( J)=R2 
GO TO 148 

B( J) = (X1*(R1 + R3)+X2*(R1 + R2) )/( ( X 1+X2 ) *2 . ) 
I F(B(J)*HCJ).LE.O.) B(J)=-.1*B(J)+EPS*R3 



C UPDATE SEARCH VARIABLES 

SI ZE a S I ZE* ( F LOAT (L)*2.)/4. 
KOUNT=KOUNT +1 
I F(NCYC.GT.N) GO TO 5 
GO TO 15 

150 SCASV 3 SCALE 

SCALE 3 SCALE/( 1 . +STEP) 

K= K+l 

SIZE** SI ZE/( l.+STEP) 

GO TO 100 

160 SIZE= S I ZE/( 1. +STEP) 

DO 180 J**1,N 
180 B( J)=H( J) 

GO TO 5 
END 


CGFM Program Variable List 


BETA 

FS 

FSS 

J 

JN 

K 

L 

NCYC 

R1,R2,R3 

SCALE 

SCASV 

STEP 

TSAVE 

TSQR 

XI, X2 


conjugate direction weighting 
saved F 
saved FS 
parameter index 
J + N 

indicator for step size reductions 

number step size increases within iteration 

number of iterations before restarting conjugate search 

terms in quadratic curve fit 

step size scale factor 

saved SCALE 

step size 

saved TSQR 

squared gradient terms sum 
partial product 
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APPENDIX C 


SAMPLE PROBLEM LISTING 


The computer listing for the optimization of a six-parameter model is presented: 


1 NAM3=( W, AMP, PHA, ND, KPR) 

2 NAM2 = ( B, I D, Kl, L 1 MIT, N, JD) 

3 I D“( 1=Z;2=G; 3*»P; 4, 5=CZD, CZW; 6, 7=CPD, CPW; 8 = DT) 

4 INST=(1=SEARCH,2=NAM2, 3=NAM3, 4= PR I NT, 5=PLANT =MODE L, 6 = ST0P) 

5 I ER=(0=CONV, l=NOT CONV, 2=ERR0R) 

6 NAM3? 


7 &nam3 w=l, 3, 7, 10, 15, 20, 25, 30, 35, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 

8 amp=l, .95, .77, .7,. 6 7,. 63,. 6,. 53, .48, .44, .3 5, .31, .33, .35, 2*. 32, . 3, .29, 

9 .27, .26, pha=-2, - 13 , -24, - 31 , -35, -44, -5 7, -62, -71, -75, -87, - 110 , -92, -10 5, 

10 -119,-128, -145, -156, -166, -172, nd=20, kpr=0 Send 

11 NAM 2? 

12 &nam2 b=l, 1, 10, 1000, 1, 1000, i d-1, 3, 3, 3, 6 , 7, kl=0, 1 i mi t = 200, n=6, jd=0 Send 

13 L IMI T, N, Kl, JD-200 6 0 0 I D= 1 3 3 3 6 7 

14 B= 1.0000E 00 1.0000E 00 1.0000E 01 1.0000E 03 1.0000E 00 1.0000E 03 

15 I NS T? 

16 1 

17 F« 5. 8507025E-01 I ER» 0 KOUNT, KNT, S I Z E= 93 417 2.868E-06 

19 ?NST? 9144E °° 3 * 6998E °° 2 * 6970E 01 3.0934E 02 2.6459E-01 1.4546E 02 

20 6 

CHCRW400 TERMINATED: STOP 


The listing shows the plant data entered in NAM3 namelist, a model built in NAM2 
namelist, and a search for the optimum parameters. The line numbers were added for 
discussion purposes, and program responses are capital letter while user input are 
lowercase. A line by line discussion of the listing follows: 


Lines 1 to 5 - For convenience the program listing of the NAM3 and NAM2 variables 
and the ID, INST, and IER variable codes 

Line 6 - A program prompt for NAM3 namelist data 

Lines 7 to 10 - The user entered system NAM3 frequency response data 


(The frequencies of W are entered in hertz rather than in radians. This requires 
associating a conversion factor of 2n with certain parameters of B. ) 


Line 11 - A program prompt for NAM2 namelist data 

Line 12 - User entered NAM2 namelist data required to build the iterative six- 
parameter transfer function model (See eq. (11). ) 

Lines 13 and 14 - Program verification of the NAM2 namelist data 

Line 15 - A program prompt for a branching instruction 
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Line 16 - A user input 1 directing the program to search for the optimum parameters 


Lines 17 and 18 - Program output of search results (Upon convergence the cost func- 
tion was reduced to 0. 585 after 93 iterations and 417 cost function evaluations. The 
optimized parameters are listed in eq. (12). ) 

Line 19 - Another program prompt for a branching instruction 

Line 20 - A user input 6 to stop the program 

In the sample problem parameter optimization the KPR variable was not set to 1; 
therefore, there was no recording of the plant data. The INST variable was not set to 
4; therefore, there was no recording of the model frequency response. The following 
listing shows the output upon recalling the program for the case KPR = 1 and INST = 4, 
but with the NAM3 and NAM2 namelist data otherwise unchanged. 


Sample problem plant and model frequency response 


NAM3=(W,AMP, PIIA / ND, KPP1 
NAM2>(B, ID, Kl, LIMIT, N,jl)) 

ID=(l=Z;2=f!;3 = P;4 / 5=CZn / CZW;B / 7=CPf\rPW; 8 = DT ) 

INST=( 1 = SEARCH, 2 = NAM2 , 3-NAM3, 4 = PR I’JT, 5 = PLA?!T=MnoF! , 6=ST0P) 
IER=(0=CONV, l=NOT CONV, 2=ERR0R ) 

NAM 3? 

Snam3 kpr=l Send 
PLANT 


w 


AMP 

PM/' 


1.00000E 

00 

1.00000E 00 

- 2 . onnnnF 

00 

3.00000E 

00 

9 . S0000E-01 

-1 , ™or>OF 

01 

7 . 00000E 

00 

7. 70000E-01 

-2. 4nnonF 

01 

1.00000E 

01 

7.00000E-01 

innonp 

01 

1.50000E 

01 

B.70000E-01 

-3 . snonnr 

01 

2.00000E 

01 

B. 30000E-01 

-4.4ononF 

01 

2.50000E 

01 

B.O00O0E-01 

-5. 7nonor 

01 

3.00000E 

01 

5. 30000E-01 

- Ft . 200H0F 

01 

3.50000E 

01 

4. 80000E-01 

-7. 1P000F 

01 

4.00000E 

01 

4.40000E-01 

-7. snonoF 

01 

5.00000E 

01 

3. 50000E-01 

- 9 . 70000E 

01 

6.00000E 

01 

3.1O000E-01 

-l. mono* 

02 

7.00000E 

01 

3.30000^-01 

20000F 

01 

8.00000E 

01 

3 . 50000E-01 

-1.H50 *0F 

02 

9.00000E 

01 

3. 20000E-01 

-1. 1 0 0 OOF 

02 

1.00000E 

02 

3 . 20000E-01 

-T . ?pnnnF 

02 

1.10000E 

0 2 

3.00000F-ni 

-1.4S0O0F 

02 

1.20000E 

02 

2.90000F-01 

-1. 5R000F 

02 

1.30000E 

02 

2. 70000E-01 

-i .nsonoF 

02 

1.40000E 

02 

2 . 60000E-01 

-1. 72000E 

02 


NAM 2? 

Snam2 Send 

LIMIT, N,Kl,JD-200 BOO JP= 1 3 3 3 F 7 

B= 4.9144E 00 3.6998E 00 2.R970E 01 3.0934E 02 2.B459F-01 1.4540E 0? 
INST? 

4 
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MODEL 


W AMP PH A 

1.00000E 00 9 . 84501E-01 -R.14033F 00 

3.00000E 00 9 . 047 3RE-01 -1.5J.B38E 01 

7.00000E 00 7.88533E-01 -2.45215E 01 

1.00000E 01 7.40270F-01 -3.01572F 01 

1.50000E 01 6.77591E-01 -3.9299RE 01 

2.00000E 01 6. 21103E-01 -4.78250F 01 

2.50000E 01 5 . 69 20 8E-01 -5.55071F 01 

3.00000E 01 5 . 22819E-01 -R.2350RF 01 

3.50000E 01 4. 82319E-01 -G.84Q4RF 01 

4.00000E 01 4 . 4749 2 E-01 -7.40447F 01 

5.00000E 01 3. 926R3E-01 -8.387R4E 01 

6 . OOOOOE 01 3. 53632F-01 -9.2RR46F 01 

7.00000E 01 3 . 26457E-01 -1.01007E 02 

8. OOOOOE 01 3 . 0832RE-01 -1.003R7E 02 

9. OOOOOE 01 2 . 97201E-01 -1.18260F 02 

1. OOOOOE 02 2.01303E-01 -1.28117E 02 

1.10000E 02 2.88359E-01 -1.39502E 02 

1.20000E 02 2.84628E-01 -1.52927E 02 

1.30000E 02 2. 74384F-01 -1.R85R2E 02 

1.40000E 02 2. 52175E-01 1.74280E 02 

F= 5. 8507025E-01 IER= 1 KOUNT, KNT,S I ZF= 0 J 1.000E-01 
B® 4.9144E 00 3.6998E 00 2.6970E 01 3.0O34E 02 2.6459E-01 7.4546E 02 
INST? 
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APPENDIX D 


TRAPEZOIDAL RULE 


The trapezoidal rule for the approximation of an integral is 


/ 


1+1 


X» — x. 

f(x) dx « — [f(x i+1 ) + f(x.)] 


This leads naturally to the series expression for a piecewise linear summation of 
trapezoidal areas 



- - x i> 


where the notation f(x.) is shorted to f. . Note that for each value of i the series re- 
quires the two values f. and f. +1 . This is reduced to a single value f. as described 
next. First, the series is written in expanded form 


2 (f 2 x 2 + f l x 2 - f 2 x l 


fjXj) + — (f^ + f 9 X^ - f 9 X 9 “ f 9 x 9 ) + . 


33 23 “ 32 " * 22 ' 


2 ( f V N d + V lXN d ' W 1 ' V 1 *^- 1 ) 


The expansion after cancellation of terms from adjacent lines is 
2 (0 + f l x 2 ' f 2 x l “ f l x l) + ~ (° + f 2 x 3 " f 3 x 2 " • • 


2 (V N d + V lXN d ' V N d- 1 ' °) 


A regrouping of the expansion forms the series 
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^l( x 2 x j) + ^2 ( x 3 ” x i) + ^3 ( x 4 " x 2^ + 


+ f N-l x N. 


- X 


N d-' 


The series is expressed in the desired form by defining x n = x, and x. T . = x 
Thus, 0 1 N d +1 N d’ 


N d 

' E f(x t )(x i+i - x i-i> 

z i=l 


Identity 

If E is assumed to be a complex function of two real variables c and d, 

E = c + jd 


it readily follows that 

(|e| 2 )’ = [(c + jd)(c - jd)]' = (c 2 + d 2 )' = 2(cc' + dd') 

For another starting point, namely, 2 He (E(E’)* ), the result is 

2 He (E (E ')* ) = 2 He [(c + jd)(c’ - jd’)] = 2(cc’ + dd’) 

Thus, the two expressions are the same, and 

(|e | 2 )’ =2 He (E (E ' )* ) 

Gradient 

The cost function gradient (eq. (6)) contains the transfer function gradient. The 
transfer function gradient in equation (7) is the transfer function times a vector Z(s, b). 

VG(s,b) = [G(s,b)][z(s,b)] (7) 

The transfer function G(s,b) is the product of factored parameter factors g m (s,b) 

(eq. (3a) to (3g)). The elements of Z(s,b) are the partial derivatives divided by the 
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corresponding factor. For example, the Z(s,b) element for a bl k factor is 



There are eight structural forms for the Z(s,b) elements corresponding to the eight 
possible parameter forms. The terms are simple, and there is only a change in sign 
between corresponding numerator and denominator factor forms. 


bl 


k : 



b2: 


b2 



b4 


: k‘ 


2s 



b5 k 

2s(b4 k ) 
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b8: -s 
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TABLE I. - SEARCH RESULTS OF DIFFERENT MODEL STRUCTURES FIT TO SAME SYSTEM DATA 


















Figure 2. - Flow chart for program MODEL. 
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